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We propose a kind of Bak-Sneppen dynamics as a general optimization technique to treat 
magnetic systems. The resulting dynamics shows self-organized criticality with power law scaling 
of the spatial and temporal correlations. An alternative method of the extremal optimization is 
also analyzed here. We provided a numerical confirmation that, for any possible value of its free 
parameter r, the extremal optimization dynamics exhibits a non-critical behavior with an infinite 
spatial range and exponential decay of the avalanches. Using the chiral clock model as our test 
system, we compare the efficiency of the two dynamics with regard to their abilities to find the 
system's ground state. 
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I. INTRODUCTION 

Nature which is composed of dead and living things is 
a system far out of equilibrium. This can be evidenced by 
measuring nature's evolution with the appropriate scales 
of space and time. In particular, the dynamics of the liv- 
ing things is basically governed by the theory of natural 
selection first proposed by Charles Darwin hi. This the- 
ory is certainly one of greatest scientific achievements of 
the nineteenth century. Very briefly, this theory simply 
says that evolution occurs not by breeding species best 
adapted to their environment but by driven to extinc- 
tion those less or poorly adapted. This mechanism leads 
to the emergence of highly specialized structures. If we 
also consider the astonishing variability of the species, we 
then can say that nature is a complex system. Indeed, for 
all wc know, nature operates at the self-organized critical 
state 0. 

One of the most fundamental characteristics of a sys- 
tem with self-organized criticality (SOC) is to exhibit a 
stationary state with a long-range power law decay of the 
spatial and temporal correlations 3j . Power law is a very 
abundant behavior found either in natural phenomena 
such as the light emitted from quasars, the earthquake's 
intensity, the water level of the Nile or as a direct result 
of human activities like the distribution of cities by size, 
the repetition of words in the Bible and in traffic jams. 

Self-organized critical systems evolve to the complex 
critical state without the interference of any external 
agent. Differently from what happens in the usual critical 
phase transition, in SOC there is not a tuning parame- 
ter. The prototypical example of SOC is a pile of sand 
Another common property of SOC is that the self- 
organization state only takes place after a very long tran- 
sient period. Last but not least, a minor change in the 
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system can cause colossal instabilities called avalanches. 
Intermittent bursts of activity separating long periods of 
quiescence is called punctuated equilibrium. Gould has 
conjectured that the biological evolution itself shows, in 
fact, this kind of equilibrium Q. 

Now, if nature can really be found in a critical self- 
organized state, is this state optimal or extremal in some 
sense?. Through the systematic rejection of the worst, 
the Darwinian evolutionary theory is a born optimiza- 
tion structure. By adding interactions and competition 
between the species we will witness the emergence of co- 
evolution. 

One model specially tailored to represent the co- 
evolutionary activities of the species is the Bak-Sneppen 
model p|. In this model, each species is located on the 
sites of a lattice and has associated a fitness value be- 
tween and 1 (randomly drawn from an uniform distri- 
bution). At each time step, the species with the smallest 
associated value as well as its nearest neighbors are se- 
lected to replace their fitness with new random numbers. 
In one dimension, after a long transient time, almost all 
species have fitness larger than the critical value 0.67 

Recently, inspired by natural processes, some heuris- 
tic optimization techniques have been proposed: genetic 
algorithms , simulated annealing and extremal op- 
timization jSi]. The extremal optimization (EO) method 
seems to be the most efficient of them since it brings 
the system faster and closer to its ground-state 0. In 
a few words, this method consists of the following rules: 
i) a discrete dynamical variable Si (initially chosen at 
random) is associated to each site i of a lattice with N 
points; ii) a fitness Xi is attributed to that site using some 
given prescription; iii) all lattice sites are then increas- 
ingly ranked according to their fitness (the site with the 
worst fitness is of rank 1); iv) a site of rank k (1 < k < N) 
is selected with probability P(k) oc k~ T (r is an arbitrary 
real positive number) and its dynamical variable Si is 
changed to S 1 ,; (S i ^ Si, compulsorily) ; v) repeat at step 
iii) as long as desired. 

In this paper, a kind of Bak-Sneppen dynamics is ap- 
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plied to a spin system which possesses a complex ground 
state structure in one dimension. The discrete fitness 
variability is identified as being responsible for the ab- 
sence of a critical self-organized state. By introducing a 
noise in the spin configuration space, we show that the 
system can now reach a critical self-organized state with 
power law correlations. This simple trick can easily be ex- 
tended to other types of discrete systems. On the other 
hand, through an analysis of the spatial and temporal 
correlations, we provided a numerical confirmation that 
the extremal optimization method is not a SOC dynam- 
ics but a non-critical algorithm with an infinite spatial 
range and an exponential time decay of the avalanches. 
The two dynamics: the Bak-Sneppen with noise (BSN) 
and the extremal optimization (EO) are then investigated 
with respect to their efficiency to minimize energy. There 
is an optimal interval of the parameter r where EO is su- 
perior to BSN. 



II. THE P-STATES CHIRAL CLOCK MODEL 
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FIG. 1: The magnetic field h versus the chirality A for the 
ground state diagram of the p — 3 chiral clock model. The 
regions < 1 >, < 12 — 13 > and < 123 > have periods 1, 
2 and 3, respectively. The marked points /, II, III and IV 
were tested for the universality of the BSN dynamics. 



To explain our main ideas, we have chosen the one di- 
mensional p-states chiral clock model (CC P ) 0, 0] as 
our experimental system. In higher dimensions, this sys- 
tem exhibits a complex phase transition diagram with 
commensurate and incommensurate phases. In one di- 
mension, the competition between the applied magnetic 
field, which tries to align the spins, and the chirality, 
which tries to flip them, gives rise to a rich ground-state 
diagram in the space of the parameters. The hamiltonian 
is given by 
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where p is a positive integer number, Si = X,..,p is the 
spin variable at the site i (with periodic boundary con- 
ditions applied) and h and A are the magnetic field and 
chirality, respectively. By symmetry arguments, h and A 
may be restricted to the intervals [0, oo] and [0, 1/2]. 

For p = 3, the chiral clock model has 3 regions with 
different ground-states (see Fig. I). We denote by < 1 >, 
< 12 — 13 > and < 123 > the regions with periods 1, 2 
and 3, respectively. The numbers correspond directly to 
the sequence of the spin states. The region < 12 — 13 > 
is super degenerate, i. e., the spin state sequence 12 
may be equally followed by 12 or 13. For a finite ring, 
this means that the degenerescence of this region grows 
exponentially with the lattice size. 



III. SEARCHING THE GROUND-STATE WITH 
A BAK-SNEPPEN DYNAMICS 

We choose the following form for the fitness Xi 
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There is a certain art and arbitrariness to select the 
fitness form. For instance, one could equally well have 
chosen a pre- factor 1/4 for the magnetic fields at the sites 
i + 1 and i — 1 and 1/2 for the site i. Different fitness 
forms may lead to quite different values of some physical 
quantities determined by the dynamics, e. g., the average 
energy in the stationary state. The fitness defined by 
Equation (2) is the form with the highest symmetry - 
the magnetic field is equally distributed between the sites 
i — 1, i and i + 1 as well as is the chirality between the 
pair of sites — 1] and [i + It represents the best 
result we got for the average energy in a bunch of trials. 

Let us now describe our procedure. Initially, to each 
site i of the ring, we attribute by chance one of the p pos- 
sible values of the spin variable Si- Using equation (2), 
we calculate all the corresponding fitness Aj and find the 
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smaller (worst) one \j. New spin variables are then ran- 
domly assigned to the sites j, j + 1 and j — So although 
the dynamics involves 3 spins, due to our fitness defini- 
tion, it affects (changes) 5 consecutive fitness (from the 
site j — 2 to j + 2). Observe that here the Bak-Sneppen 
dynamics is being applied to the spin configuration space 
not to the fitness space itself (as it was originally done 
in the evolutionary context |j|). In the final step, a new 
site with the smallest fitness is searched and the whole 
process continues as long as we wish. However, even for 
a modest lattice size, this dynamics is hampered by the 
discrete fitness variability. Let us explain why. If p = 3 
(h and A fixed) there are at most 27 possible values of the 
fitness (besides, some of these may be degenerate). This 
means that, for a lattice with a reasonable size, we will 
generally find not just one but a large number of sites 
with exactly the same (smallest) fitness value. The sim- 
plest solution to the problem seems to put all those sites 
in a list and to choose one of them at random. We will 
call this dynamics as the Bak-Sneppen with draw (BSD). 
As we shall see later, neither from the point of view of 
the spatial correlations nor from the time correlations is 
this dynamics a true SOC. 



IV. THE NOISE 



total) as the transient time, averages were then taken 
over the remaining steps. The results for the energy his- 
tograms are shown in Figure 2. 
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FIG. 2: Energy per spin histograms for the EO (with r = 2.1 
and 1.2), BSN and BSD dynamics. They were determined at 
the point III of the Figure 1. The exact ground state energy 
per spin is 0.8055. 



It is the discrete variability of the fitness (or, equiv- 
alently, of the spin variable) which precludes the estab- 
lishment of a self-organized critical state, so we searched 
a simple way to solve the problem. Suppose that a noise 
r\ is added to the original spin variable Si (an integer in 
the interval [l,f>]), i- e. 

S' i =S i + V (l-2r) (3) 

where r 6 [0, 1] is a random number generated from a 
flat distribution. The new spin variable as well as the 
corresponding fitness have now the desired continuous 
characteristic. By choosing the noise r\ sufficiently small 
the relevant physical properties (like the energy per site 
or the magnetization per site) will not be affected. In 
this paper we set r/ = 10~ 3 . Using double numerical 
precision, we varied rj from 10~ 12 up to 10~ 3 with no 
important difference for the physical quantities. With 
this trick, any previously discrete fitness A is turned into 
a continuous variable inside some interval AA controlled 
by the noise n, the chirality A and the magnetic field h. 
We will call this dynamics as the Bak-Sneppen with noise 
(BSN). Our next objective is to compare the 3 dynamics 
EO, BSD and BSN with respect to their efficiencies. 

In order to guarantee that the stationary state had 
already been reached and to get good averages for the 
physical quantities, a huge amount of computation was 
performed. Using the same fitness definition (equation 
(2)) for all the three algorithms BSD, BSN and EO, we 
simulated each one of them over 2 10 9 runs on a ring of 
4001 sites. Discarding the first 2 10 8 runs (10% of the 



The simulations were done at the point III of the re- 
gion < 123 > (see Figure 1) where the ground state en- 
ergy per spin is equal to 0.8055. The best performance 
was obtained by the EO algorithm with r = 2.1, fol- 
lowed by the EO with r = 1.2, BSN and BSD. Their re- 
spective mean energies were 0.828 ± 0.003, 0.969 ± 0.009, 
0.973 ± 0.004 and 0.983 ± 0.005. 



V. THE SPATIAL AND TEMPORAL 
DISTRIBUTIONS 

To study the spatial correlation, we measured the dis- 
tribution D(x), of the distance x between two subsequent 
activated sites y . Recall that for the three algorithms 
we activate a site i by changing its spin state Si and the 
corresponding fitness A^. We plotted D(x) in the Figure 
3. 

Clearly, EO is an algorithm with an infinite spatial 
range - it doesn't matter how far one site is of the other, 
the activation probability is constant. On the other 
hand, the BSN algorithm exhibits a power law depen- 
dence D(x) ~ x _319±0 02 with the exponent being com- 
patible with that of the Bak-Sneppen model (3.23 ± 0.02 
|Tl|). Moreover, we got at the points /, //, III and IV 
(see Figure 1) the same exponent value, indicating that 
the universality principle holds in the space of the param- 
eters h and A. Finally, from the Figure 3 we conclude 
that the BSD algorithm is of a mixed kind, having the 
BSN behavior for small distances and the EO for large 
distances. 
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FIG. 3: The probability distribution D(x) versus the distance 
x between two subsequent activated sites. The algorithm EO 
shows an infinite spatial range while the BSN algorithm has 
a power law scaling. 



In the stationary regime, the BSN dynamics shows an 
intermittent time evolution with long periods of passivity 
interrupted by sudden bursts of activity, i. e., it exhibits 
a punctuated equilibrium. This abrupt change of activ- 
ity is called an avalanche. We say that an avalanche is 
happening if there exist one site of the ring whose fitness 
is smaller than a certain threshold A c . The size A of an 
avalanche is defined as the number of subsequent time 
steps with at least one fitness below that threshold. For 
each one of the three dynamics we calculated the prob- 
ability distribution of the avalanches P(A). The results 
are illustrated in the Figure 4. 

The Figures 4(a), 4(b) and 4(c) correspond to the BSN, 
EO and BSD algorithms, respectively. Before we analyze 
P(A), let us first discuss how the critical thresholds were 
obtained. At the time step 2 10 8 (the end of the tran- 
sient period) and with the magnetic field and chirality 
fixed at the values 0.8 and 0.35, respectively, we deter- 
mined the fitness histograms of the whole ring. They are 
shown as insets of the Figure 4. Since in the EO and the 
BSD algorithms there are only some few possible values 
for the fitness A (a maximum of 27 when p — 3 and many 
of them are degenerate), we may vary the (discrete) A 
values until the avalanches are founded. Remember that 
if A < A c there will be no avalanches while if A > A c only 
one (infinite) avalanche is present. Due to the discrete- 
ness values of A, it is very improbable to find the system 
in a sub or supercritical regime. We found A c = —0.8055 
and A c = —1.5061 for the EO and BSD dynamics, re- 
spectively. For the BSN algorithm, however, the fitness 
is continuous. There is a fitness fluctuation AA due to the 
presence of the noise rj. In our case, the maximum value 
of the relative fluctuation AA / A is 1.86 10" 3 . To find 
the critical threshold, we wrote a program which varies 
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FIG. 4: Plot of the probability distribution of the avalanches 
P(A) versus the avalanches size A for the BSN, EO and BSD 
algorithms in (a), (b) and (c), respectively. The insets corre- 
spond to the fitness histograms evaluated at the time step 2 
10 s . 



the fitness values (inside the range AA) until the resulting 
avalanches have sizes between 500 and 2000 (which seems 
reasonable to avoid the sub and supercritical regimes). 
We determined A c = -1.5059. 

We can now go back to the probability distribution 
of the avalanches P(A). They were all calculated at 
their respective thresholds A c . The EO algorithm (Fig- 
ure 4(b)) shows an exponential decay with A. The 
BSN (Figure 4(a)), on the other hand, has a power law 
scaling P{A) ~ ^-0-87±o.oi Tmg va i ue j s somewhat 
discrepant from that of the evolutionary Bak-Sneppen 
model (1.07 ± 0.01 |llj), but we should remember that, 
for our system, the exact determination of the critical 
threshold position is much more difficult and one can 
easily be found in a sub or supercritical regime. For 
the BSD algorithm (Figure 4(c)), we found once more a 
mixed behavior described by P(A) ~ ^4-o.68 e -o.oo2A_ 
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VI. CONCLUSIONS 

Based on the Bak-Sneppen ideas, we proposed an opti- 
mization algorithm whose dynamics, in the steady state 
regime, exhibits self-organized criticality. It has been suc- 
cessfully applied to find the ground states of the clock 
chiral model and it can be easily extended to any other 
magnetic system. From the three dynamics studied in 
this paper, we conclude that only the BSN algorithm con- 
ducts the system to a self-organized critical state. The 
essential point was to recognize that the discrete fitness 
values, present in the EO and BSD algorithms, preclude 
the system to have a power law decay of the spatial corre- 
lation. To overcome this problem we introduced a small 
noise in the spin configuration space. Another important 
conclusion we arrived is that to possess SOC characteris- 
tics does not guarantee a dynamics to be optimal. This 
statement was proved here for the EO algorithm. The 
EO algorithm, which is not SOC since it has an infinite 
spatial range and the avalanches decay exponentially, has 
a peak of the energy histogram (with r = 2.1) which 



is only 2.7% further from the exact ground state value. 
This is a stunning result which cannot be beaten by the 
BSN dynamics. At r = 1.2, the EO and BSN are almost 
equivalent. On the other hand, the BSD algorithm is of a 
mixed type. It can be obtained by taking the limit r\ — > 
in the BSN. 

Finally, it is important to observe that, contrary to 
what happens with BSN, the EO algorithm has one ar- 
bitrary and free parameter r. If r — > 0, EO becomes a 
random walk and if r — > oo, the system may stack in a 
dead end. In both limits, the EO efficiency is completely 
spoiled and an optimal value of r should be found in 
somewhere between. In the BSN dynamics there is no 
tuning parameter. Using the ideas described here, we 
wish to develop a SOC optimization algorithm to some 
classical combinatorial problems like the bipartition of 
graphs. 
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